rm(list=ls())

# Set working directory to the repository root
# Update this path if running from a different location
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
setwd("..")
sink("EC_log.txt")



library(tidyverse)

`%notin%` = Negate(`%in%`)


respondents <- readRDS(file = "Data/EliteGWData.rds")

respondents$weight_use %>% summary()




########
print("Number of Unique Respondents")
length(unique(respondents$ResponseId))

warming <- respondents %>% dplyr::select(pid3, weight_use, climatedeny_1, climatedeny_2, climatedeny_3)

respondents <- respondents %>% mutate(climatedeny_1_bin = case_when(climatedeny_1 == "Yes" ~ 1,
                                                            is.na(climatedeny_1) ~ NA_real_,
                                                            TRUE ~ 0), 
                              climatedeny_2_bin = case_when(climatedeny_2 == "Caused mostly by human activities." ~ 1,
                                                            is.na(climatedeny_2) ~ NA_real_,
                                                            TRUE ~ 0), 
                              climatedeny_3_bin = case_when(climatedeny_3 == "Most scientists think global warming is happening." ~ 1,
                                                            is.na(climatedeny_3) ~ NA_real_,
                                                            TRUE ~ 0))


warming <- warming %>% mutate(climatedeny_1_bin = case_when(climatedeny_1 == "Yes" ~ 1,
                                                            is.na(climatedeny_1) ~ NA_real_,
                                                            TRUE ~ 0), 
                              climatedeny_2_bin = case_when(climatedeny_2 == "Caused mostly by human activities." ~ 1,
                                                            is.na(climatedeny_2) ~ NA_real_,
                                                            TRUE ~ 0), 
                              climatedeny_3_bin = case_when(climatedeny_3 == "Most scientists think global warming is happening." ~ 1,
                                                            is.na(climatedeny_3) ~ NA_real_,
                                                            TRUE ~ 0))
library(weights)
print("Weighted Chi Square Tests by party for the three climate denial questions")

round(wtd.chi.sq(warming$climatedeny_1_bin, warming$pid3, weight = warming$weight_use), digits =5)
round(wtd.chi.sq(warming$climatedeny_2_bin, warming$pid3, weight = warming$weight_use), digits =5)
round(wtd.chi.sq(warming$climatedeny_3_bin, warming$pid3, weight = warming$weight_use), digits =5)


warming_long <- warming %>% dplyr::select(pid3, weight_use, climatedeny_1_bin, climatedeny_2_bin, climatedeny_3_bin) %>% pivot_longer(cols = starts_with("clim"), names_to = "question_index", values_to = "binary_trust")

warming_cells <- warming_long %>% filter(!is.na(binary_trust)) %>% group_by(pid3, question_index, binary_trust) %>% summarise(cell_mass = sum(weight_use))
warming_margins <- warming_cells %>% group_by(pid3, question_index) %>% summarise(margin_mass = sum(cell_mass))
warming_cells <- warming_cells %>% left_join(warming_margins) %>% mutate(proportion = cell_mass/margin_mass)

warming_cells <- warming_cells %>% mutate(question_long = case_when(question_index == "climatedeny_1_bin" ~ "Do you think that global warming is happening?\n[Yes]",
                                                                    question_index == "climatedeny_2_bin" ~ "Assuming global warming is happening, do you think it is\n[caused mostly by human activities?] ",
                                                                    question_index == "climatedeny_3_bin" ~ "Which comes closest to your own view?\n[Most scientists think global warming is happening.]",
))
warming_cells$pid3 <- fct_rev(warming_cells$pid3)
warming_cells$question_long <- fct_relevel(warming_cells$question_long, "Do you think that global warming is happening?\n[Yes]", 
                                           "Assuming global warming is happening, do you think it is\n[caused mostly by human activities?] ", 
                                           "Which comes closest to your own view?\n[Most scientists think global warming is happening.]")

p1 <- ggplot(warming_cells, aes(x=pid3, y=proportion, fill=as.factor(binary_trust))) + geom_bar(stat="identity", color="white") + facet_wrap(~question_long, ncol=1) + coord_flip() +
  theme_minimal() + scale_fill_manual("Did the respondent select the scientifically\naccurate response [shown in brackets]", values = c("#7b3294", "#008837"), labels = c("No", "Yes"),  guide = guide_legend(reverse = TRUE)) +
  theme(legend.position = "top", strip.text.x = element_text(size = 12,face = "bold", hjust=0), axis.text.y = element_text(color = c( "#b2182b","#666666",  "#2166ac"), face= "bold", size=12)) +
  scale_y_continuous(labels = scales::percent) + xlab("") + ylab("Proportion of Political Elites and Public Servants")
ggsave("Fig1A.pdf",p1, width = 4, height = 6)

respondents_sci_conspiracy <- respondents %>% dplyr::select(pid3, weight_use, conspiracy_4)
respondents_sci_conspiracy_long <- respondents_sci_conspiracy %>% pivot_longer(cols = starts_with("cons"), names_to = "statement_index", values_to = "percent_likely" )

respondents_sci_conspiracy_long$pid3 <- fct_relevel(respondents_sci_conspiracy_long$pid3, "Rep", "Ind/DK", "Dem") 

library(ggridges)

respondents_sci_conspiracy_long <- respondents_sci_conspiracy_long %>% mutate(statement_long = case_when(statement_index == "conspiracy_4" ~ "The scientific consensus on global warming is based on\nmanipulated data to suppress dissent and undermine\nthe economic dominance of the United States.", 
                                                                                                         statement_index == "conspiracy_5" ~ "COVID-19 was engineered and released by a lab in Wuhan,\nChina where the pandemic originated."))

(p2 <- ggplot(respondents_sci_conspiracy_long, aes(x=percent_likely, fill = pid3, group=pid3, y=pid3)) + geom_density_ridges(aes(height=..scaled.., weight=weight_use), scale= 0.95, stat="density")+ theme_minimal() +
    facet_wrap(~statement_long, ncol = 1) + scale_fill_manual("", values=c( "#b2182b","#666666",  "#2166ac")) + theme(legend.position = "none") + xlab( "Respondent's rating of how likely it is on a [0-100] scale that the statement is true") + ylab("") +
    theme(strip.text.x = element_text(size = 12,face = "bold", hjust=0), axis.text.y = element_text(color = c( "#b2182b","#666666",  "#2166ac"), face= "bold", size=12))
)




library(sjstats)
two_party <- respondents %>% filter(pid3 %in% c("Rep", "Dem"))

mann_whitney_test(two_party, "conspiracy_4", by = "pid3", weights = "weight_use")


library(spatstat)

print("Weighted t-test of climate conspiracy by party, indepoendents excluded")
weights::wtd.t.test(x=respondents$conspiracy_4[respondents$pid3 == "Dem"], y = respondents$conspiracy_4[respondents$pid3 == "Rep"], weight =respondents$weight_use[respondents$pid3 == "Dem"], weighty =respondents$weight_use[respondents$pid3 == "Rep"], bootse = T, bootp=T, bootn=10000)


ec1 <- ewcdf(two_party$conspiracy_4[two_party$pid3 == "Rep"], two_party$weight_use[two_party$pid3 == "Rep"], normalise = T)


ec2 <- ewcdf(two_party$conspiracy_4[two_party$pid3 == "Dem"], two_party$weight_use[two_party$pid3 == "Dem"], normalise = T)

#plot(ec1, verticals=TRUE, do.points=FALSE, col = "red")
#plot(ec2, verticals=TRUE, do.points=FALSE, add=TRUE, col='blue')

print("Rep percentage less than 10 on conspiracy ECDF")
ec1(9.99999)

print("Dem percentage less than 10 on conspiracy ECDF")
ec2(9.99999)

print("Rep percentage greater than 50 on conspiracy ECDF")
1-ec1(49.99999)

print("Dem percentage greater than 50 on conspiracy ECDF")
1-ec2(49.99999)

print("Rep percentage less than 75 on conspiracy ECDF")
ec1(75)


ec_plot1 <- ggplot(data.frame(x=seq(0,100,1)), aes(x=x)) + stat_function(fun=ec1, color="red", geom= "step") + stat_function(fun=ec2, color="blue", geom= "step") + 
  theme_minimal() +xlab("Rating") + ylab("Share of Respondents\nRating Statement Below X") + scale_y_continuous(labels = scales::percent) +
   theme(text = element_text(size = 8))

library(cowplot)
fplot2 <- cowplot::plot_grid(p2, ec_plot1, labels = "AUTO", rel_heights  = c(4, 6), ncol=1)


ggsave("Fig1D_2.pdf",  fplot2, device = "pdf", width = 4, height = 6, units="in")



#### Regressions
respondents$pid <- fct_relevel(respondents$pid,"Independent", "Strong Dem", "Weak Dem", "Dem Lean",  "Rep Lean", "Weak Rep", "Strong Rep")

m_gw <- lm(AGW_yes ~pid + stdz(trust_science_mean, weight_use) + stdz(science_institutions_trust_mean, weight_use) + stdz(ideology.irt, weight_use) + stdz(sexism_mean, weight_use) + stdz(racial_resentment_mean, weight_use) + stdz(conspiracy_mean, weight_use) + Group, weight = weight_use, data = respondents)
m_gw_logit <- glm(AGW_yes ~pid + stdz(trust_science_mean, weight_use) + stdz(science_institutions_trust_mean, weight_use) + stdz(ideology.irt, weight_use) + stdz(sexism_mean, weight_use) + stdz(racial_resentment_mean, weight_use) + stdz(conspiracy_mean, weight_use) + Group, family = binomial(link = "logit"), weight = weight_use, data = respondents)

m_gw_r <- lm(AGW_yes ~stdz(trust_science_mean, weight_use) + stdz(science_institutions_trust_mean, weight_use) + stdz(ideology.irt, weight_use) + stdz(sexism_mean, weight_use) + stdz(racial_resentment_mean, weight_use) + stdz(conspiracy_mean, weight_use) + Group, weight = weight_use, data = filter(respondents, pid3 == "Rep"))
m_gw_d <- lm(AGW_yes ~stdz(trust_science_mean, weight_use) + stdz(science_institutions_trust_mean, weight_use) + stdz(ideology.irt, weight_use) + stdz(sexism_mean, weight_use) + stdz(racial_resentment_mean, weight_use) + stdz(conspiracy_mean, weight_use) + Group, weight = weight_use, data = filter(respondents, pid3 == "Dem"))

m_gw_r_logit <- glm(AGW_yes ~stdz(trust_science_mean, weight_use) + stdz(science_institutions_trust_mean, weight_use) + stdz(ideology.irt, weight_use) + stdz(sexism_mean, weight_use) + stdz(racial_resentment_mean, weight_use) + stdz(conspiracy_mean, weight_use) + Group, weight = weight_use, data = filter(respondents, pid3 == "Rep"), family = binomial(link = "logit"))
m_gw_d_logit <- glm(AGW_yes ~stdz(trust_science_mean, weight_use) + stdz(science_institutions_trust_mean, weight_use) + stdz(ideology.irt, weight_use) + stdz(sexism_mean, weight_use) + stdz(racial_resentment_mean, weight_use) + stdz(conspiracy_mean, weight_use) + Group, weight = weight_use, data = filter(respondents, pid3 == "Dem"), family = binomial(link = "logit"))

library(modelsummary)

modelsummary(list("LPM" = m_gw,  "LPM - Rep" = m_gw_r, "LPM - DEM" = m_gw_d, "Logit GLM" = m_gw_logit,  "Logit GLM - Rep" = m_gw_r_logit, "Logit GML - DEM" = m_gw_d_logit), "AGW_Models.docx")


car::vif(m_gw)
summary(m_gw)
library(sjPlot)
AGW_fullmodelplot <- plot_model(m_gw,  robust = TRUE, 
                                rm.terms = c("GroupCorporate", "GroupCourts", "GroupFederalGov", "GroupLawLobbyists", "GroupMedia", "GroupNonProfit", "GroupStateLocal"), 
                                show.values = T, show.p = T, value.offset = .3,
                                prefix.labels = "none",
                                axis.labels = c("Conspiracy Belief", "Racial Resentment", "Hostile Sexism", "Ideology", "Trust in Scientific Institutions", "General Trust in Science", "Strong Rep", "Weak Rep", "Lean Rep", "Lean Dem", "Weak Dem", "Strong Dem"),
                                group.terms = c(1,1,1,2,2,2,3,3,3,3,3,3), 
                                colors = c( "#4393c3", "#d6604d",  "#666666")) + theme_minimal() + ylim(-.50,.50) +geom_hline(yintercept = 0, linetype = 3) + ggtitle("DV: [Yes] to Global Warming is 'caused mostly by human activities.'")

af <- anova(m_gw)
afss <- af$"Sum Sq"
anove_results <- cbind(af,PctExp=afss/sum(afss)*100)
anove_results <- anove_results %>% as_tibble(rownames = "Var")
library(modelsummary)
datasummary_df(anove_results, output = "AGW_Anova.docx")




AGW_partyplot <- plot_models(m_gw_r, m_gw_d, robust = TRUE, show.values = T, show.p = T, 
                             axis.labels = c("Conspiracy Belief", "Racial Resentment", "Hostile Sexism", "Ideology", "Trust in Scientific Institutions", "General Trust in Science"),
                             rm.terms = c("GroupCorporate", "GroupCourts", "GroupFederalGov", "GroupLawLobbyists", "GroupMedia", "GroupNonProfit", "GroupStateLocal")) +
  theme_minimal() + ylim(-.5,.5) +geom_hline(yintercept = 0, linetype = 3) + ggtitle("DV: [Yes] to Global Warming is 'caused mostly by human activities.'") +
  theme(legend.position = "top", legend.title = element_blank()) + scale_color_manual(values = c( "#4393c3", "#d6604d"), labels = c("Democrats", "Republicans"))

z_vars <- c(
  "trust_science_mean",
  "science_institutions_trust_mean",
  "ideology.irt",
  "sexism_mean",
  "racial_resentment_mean",
  "conspiracy_mean"
)

# 1) Standardize BEFORE splitting (using full-sample weights)
respondents_z <- respondents %>%
  mutate(across(
    all_of(z_vars),
    ~ weights::stdz(.x, weight_use),
    .names = "{.col}_z"
  ))

respondents_z$pid <- gsub(" ", "_",respondents_z$pid  )
respondents_group <- respondents_z %>% split(.$Group)

fml <- AGW_yes ~ pid +
  trust_science_mean_z +
  science_institutions_trust_mean_z +
  ideology.irt_z +
  sexism_mean_z +
  racial_resentment_mean_z +
  conspiracy_mean_z

# Run the same model within each Group (no Group FE needed post-split)
models_by_group <- respondents_group %>%
  imap(~ lm(fml, data = .x, weights = weight_use))

library(modelsummary)

AGW_subgroup_plot <- modelplot(models_by_group) + 
  scale_color_brewer(palette = "Dark2") + 
  theme(legend.position = "top") + geom_vline(xintercept = 0, linetype = 3)

ggsave("AGWSubgroup_modelplot.pdf", AGW_subgroup_plot, width = 7, height = 7)
# function to run the restricted test on a single model
pid_equality_test <- function(mod,
                              lhs = "pidStrong_Dem",
                              rhs = "pidStrong_Rep") {
  
  # make sure both coefficients exist
  coefs <- names(coef(mod))
  if (!(lhs %in% coefs && rhs %in% coefs)) {
    return(tibble(
      term = paste(lhs, "=", rhs),
      statistic = NA_real_,
      p.value = NA_real_,
      df = NA_real_
    ))
  }
  
  lh <- car::linearHypothesis(
    mod,
    paste0("", lhs, " = ", rhs, ""),
    test = "F"
  )
  
  tibble(
    term = paste(lhs, "=", rhs),
    statistic = lh$F[2],
    p.value = lh$`Pr(>F)`[2],
    df = lh$Df[2]
  )
}

# apply across model list
pid_tests <- models_by_group %>%
  imap_dfr(
    ~ pid_equality_test(.x) %>% mutate(Group = .y)
  )

write_csv(pid_tests, "SubgroupAGWRegressionStrongPartisanFtests.csv")



respondents <- respondents %>% dplyr::mutate(conspiracy_4_rev = 100 - conspiracy_4)

m_gwc <- lm(stdz(conspiracy_4_rev, weight_use) ~pid + stdz(trust_science_mean, weight_use) + stdz(science_institutions_trust_mean, weight_use) +  stdz(ideology.irt, weight_use) + stdz(sexism_mean, weight_use) + stdz(racial_resentment_mean, weight_use) + stdz(conspiracy_mean, weight_use) + Group, weight = weight_use, data = respondents)
m_gwc_r <- lm(stdz(conspiracy_4_rev, weight_use) ~ stdz(trust_science_mean, weight_use) + stdz(science_institutions_trust_mean, weight_use) + stdz(ideology.irt, weight_use) + stdz(sexism_mean, weight_use) + stdz(racial_resentment_mean, weight_use) + stdz(conspiracy_mean, weight_use) + Group, weight = weight_use, data = filter(respondents, pid3 == "Rep"))
m_gwc_d <- lm(stdz(conspiracy_4_rev, weight_use) ~ stdz(trust_science_mean, weight_use) + stdz(science_institutions_trust_mean, weight_use) +  stdz(ideology.irt, weight_use) + stdz(sexism_mean, weight_use) + stdz(racial_resentment_mean, weight_use) + stdz(conspiracy_mean, weight_use) + Group, weight = weight_use, data = filter(respondents, pid3 == "Dem"))
car::vif(m_gwc)

af <- anova(m_gwc)
afss <- af$"Sum Sq"
anove_results <- cbind(af,PctExp=afss/sum(afss)*100)
anove_results <- anove_results %>% as_tibble(rownames = "Var")
library(modelsummary)
datasummary_df(anove_results, output = "GWC_Anova.docx")

GWC_fullmodelplot <- plot_model(m_gwc,  robust = TRUE, 
                                rm.terms = c("GroupCorporate", "GroupCourts", "GroupFederalGov", "GroupLawLobbyists", "GroupMedia", "GroupNonProfit", "GroupStateLocal"), 
                                show.values = T, show.p = T, value.offset = .3,
                                prefix.labels = "none",
                                axis.labels = c("Conspiracy Belief", "Racial Resentment", "Hostile Sexism", "Ideology", "Trust in Scientific Institutions", "General Trust in Science", "Strong Rep", "Weak Rep", "Lean Rep", "Lean Dem", "Weak Dem", "Strong Dem"),
                                group.terms = c(1,1,1,2,2,2,3,3,3,3,3,3), 
                                colors = c( "#4393c3", "#d6604d",  "#666666")) + theme_minimal() +geom_hline(yintercept = 0, linetype = 3) + ggtitle("DV: 100 - Belief in Global Warming Science Conspiracy [0-100] (standardized to mu = 0, sd =1)") + ylim(-1,1)

GWC_partyplot <- plot_models(m_gwc_r, m_gwc_d, robust = TRUE,show.values = T, show.p = T, 
                             rm.terms = c("GroupCorporate", "GroupCourts", "GroupFederalGov", "GroupLawLobbyists", "GroupMedia", "GroupNonProfit", "GroupStateLocal"),
                             axis.labels = c("Conspiracy Belief", "Racial Resentment", "Hostile Sexism", "Ideology", "Trust in Scientific Institutions", "General Trust in Science")
) + theme_minimal() + ylim(-1,1) +geom_hline(yintercept = 0, linetype = 3) + ggtitle("DV: 100 - Belief in Global Warming Science Conspiracy [0-100] (standardized to mu = 0, sd =1)") +
  theme(legend.position = "top", legend.title = element_blank()) + scale_color_manual(values = c( "#4393c3", "#d6604d"), labels = c("Democrats", "Republicans"))


Reg_grid <- cowplot::plot_grid(AGW_fullmodelplot,GWC_fullmodelplot,AGW_partyplot, GWC_partyplot)

ggsave("GWRegPlotgrid.pdf", Reg_grid, device = "pdf", width = 6, height = 6)



modelsummary(models = list("Anthropogenic Global Warming" = m_gw,
                           "Anthropogenic Global Warming Party (Republicans)" = m_gw_r, 
                           "Anthropogenic Global Warming Party (Democrats)" = m_gw_d, 
                           "Global Warming Conspiracy" = m_gwc,
                           "Global Warming Conspiracy (Republican)" = m_gwc_r,
                           "Global Warming Conspiracy (Democrats)" = m_gwc_d), output = "ClimateRegModels.docx")

modelsummary(models = list("Logit Anthropogenic Global Warming" = m_gw_logit,
                           "Logit Anthropogenic Global Warming Party (Republicans)" = m_gw_r_logit, 
                           "Logit Anthropogenic Global Warming Party (Democrats)" = m_gw_d_logit), 
             output = "LogitClimateRegModels.docx")


respondents <- respondents %>% mutate(climatedeny_1_bin = case_when(climatedeny_1 == "Yes" ~ 1,
                                                            is.na(climatedeny_1) ~ NA_real_,
                                                            TRUE ~ 0), 
                              climatedeny_2_bin = case_when(climatedeny_2 == "Caused mostly by human activities." ~ 1,
                                                            is.na(climatedeny_2) ~ NA_real_,
                                                            TRUE ~ 0), 
                              climatedeny_3_bin = case_when(climatedeny_3 == "Most scientists think global warming is happening." ~ 1,
                                                            is.na(climatedeny_3) ~ NA_real_,
                                                            TRUE ~ 0))


library(fixest)
respondents
#m_gwc <- feols(stdz(conspiracy_4, weight_use) ~ climatedeny_1_bin + climatedeny_2_bin + climatedeny_3_bin + pid + stdz(trust_science_mean, weight_use) + stdz(science_institutions_trust_mean, weight_use) +  stdz(ideology.irt, weight_use) + stdz(sexism_mean, weight_use) + stdz(racial_resentment_mean, weight_use) + stdz(conspiracy_mean, weight_use) | Group, weights = ~weight_use, data = respondents)
#m_gwc_r <- feols(stdz(conspiracy_4, weight_use) ~ climatedeny_1_bin + climatedeny_2_bin + climatedeny_3_bin + pid + stdz(trust_science_mean, weight_use) + stdz(science_institutions_trust_mean, weight_use) +  stdz(ideology.irt, weight_use) + stdz(sexism_mean, weight_use) + stdz(racial_resentment_mean, weight_use) + stdz(conspiracy_mean, weight_use) | Group, weights = ~weight_use, data  = filter(respondents, pid3 == "Rep"))
#m_gwc_d <- feols(stdz(conspiracy_4, weight_use) ~ climatedeny_1_bin + climatedeny_2_bin + climatedeny_3_bin + pid + stdz(trust_science_mean, weight_use) + stdz(science_institutions_trust_mean, weight_use) +  stdz(ideology.irt, weight_use) + stdz(sexism_mean, weight_use) + stdz(racial_resentment_mean, weight_use) + stdz(conspiracy_mean, weight_use) | Group, weights = ~weight_use, data  = filter(respondents, pid3 == "Dem"))

#library(modelsummary)

#modelsummary(list("All" = m_gwc, "Rep" = m_gwc_r, "Dem" = m_gwc_d), output = "FixestGWCRegModels.docx")

sink()
respondents_z$pid <- fct_relevel(respondents_z$pid,"Independent", "Strong_Dem", "Weak_Dem", "Dem_Lean",  "Rep_Lean", "Weak_Rep", "Strong_Rep")

respondents_z$climate_conspiracy_z <- stdz(respondents_z$conspiracy_4, respondents_z$weight_use)
respondents_group <- respondents_z %>% split(.$Group)

fml_clim <- climate_conspiracy_z ~ pid +
  trust_science_mean_z +
  science_institutions_trust_mean_z +
  ideology.irt_z +
  sexism_mean_z +
  racial_resentment_mean_z +
  conspiracy_mean_z

# Run the same model within each Group (no Group FE needed post-split)
models_by_group_clim <- respondents_group %>%
  imap(~ lm(fml_clim, data = .x, weights = weight_use))

library(modelsummary)

GWC_subgroup_plot <- modelplot(models_by_group_clim) + 
  scale_color_brewer(palette = "Dark2") + 
  theme(legend.position = "top") + geom_vline(xintercept = 0, linetype = 3)


ggsave("GWCSubgroup_modelplot.pdf", GWC_subgroup_plot, width = 7, height = 7)
# function to run the restricted test on a single model
pid_equality_test <- function(mod,
                              lhs = "pidStrong_Dem",
                              rhs = "pidStrong_Rep") {
  
  # make sure both coefficients exist
  coefs <- names(coef(mod))
  if (!(lhs %in% coefs && rhs %in% coefs)) {
    return(tibble(
      term = paste(lhs, "=", rhs),
      statistic = NA_real_,
      p.value = NA_real_,
      df = NA_real_
    ))
  }
  
  lh <- car::linearHypothesis(
    mod,
    paste0("", lhs, " = ", rhs, ""),
    test = "F"
  )
  
  tibble(
    term = paste(lhs, "=", rhs),
    statistic = lh$F[2],
    p.value = lh$`Pr(>F)`[2],
    df = lh$Df[2]
  )
}

# apply across model list
pid_tests_clim <- models_by_group_clim %>%
  imap_dfr(
    ~ pid_equality_test(.x) %>% mutate(Group = .y)
  )

write_csv(pid_tests, "SubgroupGWCRegressionStrongPartisanFtests.csv")



#group DK and IND together
respondents$pid3[respondents$pid3 == "DK"] <- "Ind"

#### Generate Summary Statistics of Respondents
library(gtsummary)
library(xtable)
sink("SI_Table1.tex")
xtable(questionr::wtd.table(respondents$Group, weights = respondents$weight_use))
sink()
respondents_dec <- respondents %>% dplyr::select(pid3, 
                                                 race,
                                                 hisp,
                                                 gender,
                                                 lgbt,
                                                 age,
                                                 Group,
                                                 weight_use)

library(survey)
respondents_dec$hisp <- fct_relevel(respondents_dec$hisp, "No", "Yes")
respondents_dec$party.relevance_1 <- fct_relevel(respondents_dec$party.relevance_1, "Strongly disagree", "Somewhat disagree", "Neither agree nor disagree", "Somewhat agree", "Strongly agree")
respondents_dec$party.relevance_2 <- fct_relevel(respondents_dec$party.relevance_2, "Strongly disagree", "Somewhat disagree", "Neither agree nor disagree", "Somewhat agree", "Strongly agree")
respondents_dec$party.relevance_3 <- fct_relevel(respondents_dec$party.relevance_3, "Strongly disagree", "Somewhat disagree", "Neither agree nor disagree", "Somewhat agree", "Strongly agree")
respondents_dec$age <- fct_relevel(respondents_dec$age, "21 - 25", "26 - 30", "31 - 35", "36 - 40", "41 - 45","46 - 50","51 - 55","56 - 60","> 60")

### recode democraphic variables for descriptive tables
respondents_dec$race2 <- "POC"
respondents_dec$race2[respondents_dec$race == "White"] <- "White"
respondents_dec$race2[is.na(respondents_dec$race)] <- NA

respondents_dec$gender2 <- "Trans/Non-Binary/Self Describe"
respondents_dec$gender2[respondents_dec$gender == "Man"] <- "Man"
respondents_dec$gender2[respondents_dec$gender == "Woman"] <- "Woman"
respondents_dec$gender2[respondents_dec$gender == "Prefer not to say"] <- NA
respondents_dec$gender2[is.na(respondents_dec$gender)] <- NA

respondents_dec$lgbt2 <- "No"
respondents_dec$lgbt2[respondents_dec$lgbt == "Yes"] <- "Yes"
respondents_dec$lgbt2[respondents_dec$lgbt == "Prefer not to say"] <- "Prefer not to say"
respondents_dec$lgbt2[is.na(respondents_dec$lgbt)] <- NA
respondents_dec <- respondents_dec %>% dplyr::select(-race, -gender, -lgbt, -Group)


### Output political elite respondents descriptive table as LaTeX
respondents_dec_svy<-svydesign(id = ~1, weights = ~weight_use, data = respondents_dec) 


table2SI <- respondents_dec_svy %>% tbl_svysummary() %>% as_gt() %>% gt::as_latex()
sink("SI_Table2.tex")
table2SI[1]
sink()

#manually add in Hispanic summary data because tbl_summary doesn't accept it 
questionr::wtd.table(respondents_dec$hisp, weights = respondents_dec$weight_use)

sink("SI_Table2_HispanicAddition.txt")
prop.table(questionr::wtd.table(respondents_dec$hisp, weights = respondents_dec$weight_use))
sink()

library(readxl)

pub_opp <- read_csv("Data/public_opinion_output.csv")

pub_opp %>% ggplot(aes(x=pid, y = mean, ymin = lower.ci, ymax = upper.ci, color = pid)) +
  geom_errorbar(width = 0, linewidth = 2) +geom_point(size = 3) + theme_minimal() +
  scale_color_manual(values =c( "#2166ac", "#666666" , "#b2182b") ) +
  theme(legend.position = "none", strip.text.y = element_text(size = 12,face = "bold", hjust=0), axis.text.x = element_text(color = c( "#2166ac", "#666666" , "#b2182b"), face= "bold", size=12)) +
  facet_wrap(~question, nrow = 1) + xlab("") + ylab("Predicted percent agreeing") + scale_y_continuous(limits = c(0,1), labels = scales::percent) -> pub_opp_plot

ggsave("PublicOpinionClimate.pdf", pub_opp_plot, width = 10, height = 4)
sink()


